threshold and dry weight temperatures are region dependent – how well does this translate to decades in the future when regional relationships may very well change?
This code block will create threshold and dry weight factor estimates for brickman points.
Caroline’s dynamic thresholds create cfin dry weight estimates based on June and August. Therefore, we will read these months and downsample by 75% to improve runtime.
bathy <- brickman::read_brickman(scenario = "PRESENT",
vars = c("Bathy_depth", "land_mask"),
interval = "ann",
form = "stars")
sst <- brickman::read_brickman(scenario="PRESENT",
vars = "SST",
interval = "mon",
form = "stars")
brickman <- c(bathy, slice(sst, month, 6), slice(sst, month, 8)) |>
st_downsample(n = 1) |>
st_as_sf(coords = c("x", "y"), crs = 4326) |>
filter(land_mask == 1) |>
na.omit() |>
rename(SST_jun = SST, SST_aug = SST.1)
Next, we map these brickman points to regions provided by Caroline. Points outside of the regions are omitted.
regions <- read_sf(dsn = "shapefiles/Regions_dw_vd_poly_all.shp") |>
st_make_valid() |>
st_transform(crs = 4326)
brickman <- mutate(brickman, region_tc =
st_intersects(brickman, regions, sparse = FALSE) |>
apply(1, function(u) ifelse(any(u), regions$id[which(u)], NA))) |>
na.omit()
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
Next, we implement Caroline’s dry weight conversion equation to determine C. Finmarchicus dry weight multiplier at each present-day brickman point. Dry weight in is units of micrograms (\(\mu g\)).
Original conversion from Caroline:
pres |>
mutate(Cfin_C5C =
ifelse(REGION_predictions %in% c("GoM", "WSS", "ESS") & MONTH %in% c(5,6),
(-18.8 * T0_50)+332,
ifelse(REGION_predictions %in% c("neGSL","nwGSL", "sGSL", "sNL")
& MONTH %in% c(7,8), (-18.8 * T0_50)+332, NA)),
Chyp_C4C =
ifelse(REGION_predictions %in% c("GoM", "WSS", "ESS") & MONTH %in% c(4,5),
((-4.71 * T0_50)+94)*6,
ifelse(REGION_predictions %in% c("neGSL","nwGSL", "sGSL", "sNL")
& MONTH %in% c(5,6), ((-4.71 * T0_50)+94)*6, NA)),
Cgla_C4C =
ifelse(REGION_predictions %in% c("GoM", "WSS", "ESS") & MONTH %in% c(4,5),
((-4.71 * T0_50)+94)*2,
ifelse(REGION_predictions %in% c("neGSL", "nwGSL", "sGSL", "sNL")
& MONTH %in% c(5,6),((-4.71 * T0_50)+94)*2, NA)))
Implementing for Brickman:
brickman <- brickman |>
rowwise() |>
mutate(idw = (-18.8 * ifelse(region_tc %in% c("GoM", "WSS", "ESS"),
SST_jun, SST_aug)) + 332) |>
ungroup()
quantile(brickman$idw, probs = c(0, .05, .50, .95, 1.00))
## 0% 5% 50% 95% 100%
## -58.07668 10.57657 98.94075 212.36640 231.51398
Expected thresholds from Caroline:
Expected dry weight values
Differences are likely due to the fact that I’m using SST instead of T0-50. Currently waiting on T0-50 values to consider, but in the meantime:
brickman <- brickman |>
mutate(bidw = map(idw + 120, ~max(82, .x)) |> unlist())
Next, we find the threshold values for the brickman points. Threshold depends on bathymetry, month, and region and is determined via the text tables provided by Caroline. We will determine threshold for the months of January, April, July, and November.
Some regions are clumped together for threshold conversions.
region_thresholds <- function(region_tc) {
if (region_tc == "GoM") {
"gom"
} else if (region_tc == "WSS") {
"wss"
} else if (region_tc %in% c("ESS", "sNL")) {
"ess"
} else {
"gsl"
}
}
brickman <- select(brickman, -land_mask, -SST_aug, -SST_jun) |>
mutate(Jan = 1, Apr = 4, Jul = 7, Nov = 11) |>
rowwise() |>
mutate(region_thresholds = region_thresholds(region_tc))
Not all months perform well in the text conversion due to lack of data, so we map to adjacent months using mapping provided by Caroline.
month_map <- matrix(c(
11, 11, 4, 4, 6, 6, 6, 9, 9, 10, 11, 11,
11, 11, 4, 4, 5, 6, 6, 9, 9, 10, 11, 11,
11, 11, 3, 5, 5, 6, 7, 8, 9, 11, 11, 11,
1, 1, 4, 4, 4, 7, 7, 8, 8, 12, 12, 12),
nrow = 12,
dimnames = list(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), c("ess", "wss", "gsl", "gom"))
)
brickman <- brickman |>
mutate(across(Jan:Nov, ~month_map[[.x, region_thresholds]]))
brickman_mon <- brickman
thresholds <- read.table("Calanus_threshold.txt", header = TRUE) |>
group_by(month, region_vd)
tkeys <- group_keys(thresholds)
tlist <- group_split(thresholds, .keep = FALSE) |>
map(~as.data.table(.x) |>
setkey(depth))
names(tlist) <- paste(tkeys$region_vd, tkeys$month)
get_threshold <- function(mon, bdepth, region) {
tlist[[paste(region, mon)]][data.table(bdepth), roll = "nearest"]$Min_cfin_mgm2
}
brickman <- brickman |>
mutate(across(Jan:Nov, ~get_threshold(.x, Bathy_depth, region_thresholds)))
## Warning: Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
Using new set of thresholds from Caroline, using minimum estimates.
thresholds <- read.table("Calanus_threshold_resting.txt", header = TRUE) |>
group_by(month, region_vd)
tkeys <- group_keys(thresholds)
tlist <- group_split(thresholds, .keep = FALSE) |>
map(~as.data.table(.x) |>
setkey(depth))
names(tlist) <- paste(tkeys$region_vd, tkeys$month)
get_threshold <- function(mon, bdepth, region) {
tlist[[paste(region, mon)]][data.table(bdepth), roll = "nearest"]$Min_cfin_mgm2
}
brickman_r <- brickman_mon |>
mutate(across(Jan:Nov, ~get_threshold(.x, Bathy_depth, region_thresholds)))
## Warning: Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
what is the necessary abundance based on the bandaid temperature and resting thresholds?
needed_ind_b <- brickman_r |>
mutate(needed_ind = (Jul * 1000)/bidw)
Version expressing needed abundance as percentile relative to actual
dataset.
root <- "/mnt/ecocast/projectdata/calanusclimate/src"
ae <- readr::read_csv(file.path(root, "tc_datasets/ae_regionvd_cfin.csv.gz"),
col_types = readr::cols())
quantile(ae$abundance, probs = c(0, .05, .5, .95, 1))
## 0% 5% 50% 95% 100%
## 0.000000e+00 3.219821e+01 3.148560e+03 6.040407e+04 1.338762e+06
Plot considering only regions with threshold less than absolute
maximum